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We investigate shear-induced crystallization in a very dense flow of mono-disperse inelastic hard 
spheres. We consider a steady plane Couette flow under constant pressure and neglect gravity. We 
assume that the granular density is greater than the melting point of the equilibrium phase diagram 
of elastic hard spheres. We employ a Navier-Stokes hydrodynamics with constitutive relations all of 
which (except the shear viscosity) diverge at the crystal packing density, while the shear viscosity 
diverges at a smaller density. The phase diagram of the steady flow is described by three parameters: 
an effective Mach number, a scaled energy loss parameter, and an integer number m: the number 
of half-oscillations in a mechanical analogy that appears in this problem. In a steady shear flow 
the viscous heating is balanced by energy dissipation via inelastic collisions. This balance can have 
different forms, producing either a uniform shear flow or a variety of more complicated, nonlinear 
density, velocity and temperature profiles. In particular, the model predicts a variety of multi-layer 
two-phase steady shear flows with sharp interphase boundaries. Such a flow may include a few zero- 
shear (solid-like) layers, each of which moving as a whole, separated by fluid-like regions. As we are 
dealing with a hard sphere model, the granulate is fluidized within the "solid" layers: the granular 
temperature is non-zero there, and there is energy flow through the boundaries of the "solid" layers. 
A linear stability analysis of the uniform steady shear flow is performed, and a plausible bifurcation 
diagram of the system, for a fixed m, is suggested. The problem of selection of m remains open. 

PACS numbers: 45.70.Mg, 83.50.Ax 



I. INTRODUCTION 

In spite of extensive experimental and theoretical stud- 
ies of dense granular flows, a theoretical description of 
these flows remains challenging yj. Multi-particle con- 
tacts and friction, intrinsic in slow dense flows, invali- 
date the kinetic theory @, El- Furthermore, even rapid 
dense flows (that is, flows dominated by binary collisions) 
present significant difficulties for analysis. Dilute and 
moderately dense flows of mono-disperse inelastic hard 
sphere fluids are describable, for not too high inelasticity 
of collisions, by Navier-Stokes hydrodynamics 3], which 
can be derived in a systematic way from the kinetic the- 
ory: a one-particle kinetic equation properly generalized 
to account for inelastic collisions [2, 01 ■ Hydrodynamic 
equations, that is conservation laws for the mass and mo- 
mentum of the granulate, and a balance equation for the 
energy, may still be valid, for a hard sphere fluid, at 
higher densities, after the disorder-order transition oc- 
curs. However, the respective constitutive relations are 
not derivable from a kinetic equation anymore. 

In this work we attempt to address this difficulty and 
suggest a possible hydrodynamic description of a sheared 
rapid granular How that exhibits crystallization. Shear- 
induced ordering in a dense granular medium has at- 
tracted much recent attention. It was investigated ex- 
perimentally (usually for slow flows) by many groups 
U [H El ■ The crystallization dynamics in inelastic hard 
sphere fluids has been also extensively studied in MD 
simulations 0, @ • In this work we will be dealing with 
an idealized model of inelastic particle collisions charac- 



terized by a constant coefficient of normal restitution r, 
and focus on a plane shear How. 

The present work is a next step in a series of recent at- 
tempts of extending granular hydrodynamics of inelastic 
hard sphere fluids to high densities |qL lloL . Grossman 
et al. 9] investigated a prototypical system of inelastic 
hard disks at zero gravity, placed in a two-dimensional 
rectangular box, one wall of which serving as a "ther- 
mal" wall. Grossman et al. suggested an equation of 
state, granular heat conductivity and inelastic energy loss 
rate which interpolated between the dilute limit and the 
close vicinity of the hexagonal close packing, where free- 
volume arguments are available. The density profiles, ob- 
tained by solving the hydrostatic equations numerically, 
were in good agreement with the results of RID simula- 
tions [jj. 

Meerson et al. y~l| considered a similar two- 
dimensional granular system, but with gravity. The sys- 
tem was driven from below by a "thermal" base. Employ- 
ing the constitutive relations suggested by Grossman et 
al. 0, Meerson and coworkers found steady-state den- 
sity and temperature profiles and observed good agree- 
ment with MD simulations, including the region of the 
"levitating cluster", where the density is very close to 
that of hexagonal close packing yT| . 

Bocquet et al. |1C| employed hydrodynamic equations 
to model a granular shear flow where the density ap- 
proached the random close packing density. They also 
compared their theory with experiment in a circular Cou- 
ette flow yjj. As a shear flow was present, Bocquet and 
coworkers had to specify, in addition to the rest of the 
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constitutive relations, the coefficient of shear viscosity. 
In experiment the granular temperature had been found 
to decrease more slowly, with an increase of the distance 
from the shear surface, than the velocity. To account 
for this finding, Bocquet et al. assumed that the shear 
viscosity diverges more rapidly, at random close packing, 
than the rest of the transport coefficients. 

We employ in this work Navier-Stokes granular hy- 
drodynamics for a description of a steady crystallized 
shear flow of an assembly of mono-disperse inelastic hard 
spheres under constant pressure. Our approach is sim- 
ilar to that taken by Bocquet et al. .10], with two im- 
portant differences. First, we focus on the ordered crys- 
talline phase which ends at the crystal (hexagonal or face- 
centered cubic) close packing density (ftfcci whereas Boc- 
quet et al. considered the metastable disordered phase 
which ends at the random close packing density 4> rcp . 
Second, we assume that the shear viscosity diverges at a 
smaller density than the rest of transport coefficients, see 
below. This assumption brings about the possibility of a 
two-phase steady flow. In general, a granular shear flow 
reaches a steady state when the viscous heating makes 
up for the energy dissipation via inelastic collisions. As 
we show here, this balance can be achieved in different 
ways, producing either a uniform shear flow, or a variety 
of flows with nonlinear density, velocity and temperature 
profiles. Working in the range of densities beyond the 
melting point, we determine the phase diagram of dense 
steady flows in terms of three parameters: an effective 
Mach number, a scaled energy loss parameter, and an 
integer number m: the number of half-oscillations in a 
mechanical analogy that appears in this problem. There 
are regions on this phase diagram where two or more 
different steady flow solutions are possible for the same 
values of the parameters. To get an additional insight, we 
perform a linear stability analysis of the uniform steady 
shear flow. Based on these results, we suggest a plau- 
sible bifurcation diagram of the system which, in some 
region of the parameter space, describes bistability and 
hysteresis. 

The rest of the paper is organized as follows. In Sec. 
2 we introduce our shear-induced crystallization scenario 
and the governing hydrodynamic equations and consti- 
tutive relations. Section 3 describes our model of a zero- 
gravity constant-pressure shear flow and focuses on the 
analysis of a crystallized steady shear flow in different 
regimes. Section 4 presents a linear stability analysis of 
the uniform shear flow solution and suggests a plausible 
bifurcation diagram of the system. Details of the linear 
stability analysis are presented in the Appendix. Section 
5 includes a brief discussion and summary of our results. 



II. PRESSURE-DENSITY DIAGRAM AND 
HYDRODYNAMIC EQUATIONS 

Figure n depicts, in the coordinates "volume fraction - 
pressure" , the phase diagram of a homogeneous macro- 



scopic system of elastic hard spheres The volume 

fraction (ft = {■n/Q)d 3 n (where n is the particle number 
density, and d is the particle diameter) changes from 
zero to (ftjee = \/27r/6 ~ 0.74, the density of crystal 
(either hexagonal, or face-centered cubic) close packing, 
the densest possible packing of spheres in three dimen- 
sions. The phase diagram includes four branches. The 
disordered, or gas/liquid branch starts at (ft = and con- 
tinues until the freezing point that occurs at (ft ~ 0.494. 
Then this branch splits into two branches. The horizon- 
tal (constant pressure) branch describes coexistence of 
the disordered and ordered phases. It starts at the freez- 
ing point and ends at the melting point at (ft ~ 0.545. 
At larger volume fractions the system is in the ordered 
crystalline phase that ends at the density of the crystal 
close packing (ftf cc . The last branch is a metastable ex- 
tension of the gas/fluid branch. It starts at the freezing 
point and ends at random close packing at (ft rcp ~ 0.64. 
Clearly, each branch is described by a separate equation 
of state. 

The phase diagram, presented in Fig. appears in the 
context of a homogeneous system in either equilibrium, 
or a metastable state. Granular systems are usually inho- 
mogeneous and, no less important, they are intrinsically 
far from equilibrium, due to inelastic collisions between 
particles. We will deal, however, with a small inelas- 
ticity of particle collisions and, in the spirit of kinetic 
theory 0,0, 

assume that the system is everywhere close 
to local thermodynamic equilibrium, so the phase dia- 
gram (and the constitutive relations presented below) are 
valid locally. Let the fluid density be sufficiently large, so 
that the volume fraction is everywhere above the melt- 
ing point. There are two possible phases here. One of 
them is the disordered phase (the the metastable branch 
in Fig. [TJ, the other one is the crystallized phase. Each 
of the two phases, disordered and ordered, can have in- 
homogeneous temperature and density profiles, without 
violating the constancy of the pressure. Furthermore, do- 
mains of disordered phase can in principle coexist with 
domains of the ordered phases (again, without violating 
the constancy of the pressure). Shear-induced crystal- 
lization apparently represents a (non-equilibrium) phase 
transition, so that the metastable disordered branch gives 
way, everywhere, to the stable ordered branch. 

This work addresses very dense flows, where the vol- 
ume fraction is close to that of the crystalline close pack- 
ing, <ftfcc ~ <ft "C (ftfcc- We assume that the ordered crys- 
talline branch of the phase diagram (see Fig. ^) has al- 
ready won in the competition with the disordered branch. 
To describe a flow in this phase, one needs hydrodynamic 
equations. These represent the mass conservation 



the momentum conservation 

n— =V P + ng, (2) 
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FIG. 1: Phase diagram of elastic hard sphere fluid, see text. 



and the energy balance 

3 dT „ n 
-ra— = -V • Q 
2 dt 



P : Vv - r , 



(3) 



where T is the energy loss rate due to the inelasticity of 
binary collisions. Here n(r,t) is the number density of 
grains, T(r, t) is the granular temperature, v(r, t) is the 
mean flow velocity, P is the stress tensor, g is the gravity 
acceleration, Q is the heat flux, and d/dt — d/dt + v • V 
is the total derivative. In the following we will put the 
particle mass to unity. The stress tensor P can be written 
as 

P = [-p(n, T) + n{n, T) tr (D )] I + 2r/(n, T) D , (4) 
where 

D = (1/2) [Vv + (Vvf] (5) 
is the rate of deformation tensor, 



D = D - - tr (D ) I 



(6) 



is the deviatoric part of D, and I is the identity tensor. 
We assume that the heat flux Q is given by the Fourier 
law 



Q = -k(u,T) VT. 



(7) 



Similarly to dilute 0] and moderately dense Q rapid 
granular flows, there can be an additional term in Eq. J7J, 
proportional to the density gradient. This term must 
vanish as r — > 1, and it can be neglected in the nearly 
elastic limit g< 1 that we are interested in throughout 
this paper. 

To make the formulation complete, one needs constitu- 
tive relations: the equation of state p = p(n, T) and the 
dependences of the transport coefficients /i, 77 and k, and 
of the energy loss rate T on 71 and T. For small and mod- 
erate densities, these relations can be derived from the 
Boltzmann or Enskog equation, properly generalized to 



account for inelastic collisions 0, |j| . For very large den- 
sities, that we are interested in, one can use free volume 
arguments 0]. The resulting equation of state is 



njT 



P = Pi 



(8) 



where n c = V2/d 3 is the particle number density at close 
packing, and pi is a numerical factor of order unity. Here 
and in the following we will write n c instead of n f cc (face- 
centered cubic) . The temperature dependence in Eq. © 
is exact. We assume that all transport coefficients (ex- 
cept for the shear viscosity 77, see below), diverge like 
(n c — n)^ 1 near the close packing density 0. Therefore, 
one can write the bulk viscosity coefficient /i(n,T), the 
thermal conductivity K,(n,T), and the energy loss rate 
r(n, T) in the following form 



,2*1/3 



M = Mi 



r = ri (1 



(n c — n) d 2 ' 

nc r 1/2 

(?7 C — 71 ) d? ' 

(n c — n) d 



(9) 



Equations Q are valid in the limit n c — n <C 77 c . The tem- 
perature dependences are exact. The numerical factors 
Hi, Ki, and T\ are of order unity and presently unknown; 
they can be found in MD simulations. The same type of 
divergence of the transport coefficients (again, except for 
the coefficient of shear viscosity 77) was assumed in Ref. 
|10|| , but in the vicinity of the random close packing den- 
sity. 

There is a recent evidence in the literature that the 
shear viscosity coefficient 77 of elastic hard disk fluid grows 
faster with the density, at high densities, than the rest of 
the transport coefficients [l6j]- We believe that this be- 
havior remains qualitatively correct also for three dimen- 
sional systems. Bocquet et al. accounted for this fact 
in their description of the plane shear flow near the ran- 
dom close packing density. They proposed, by analogy 
with the behavior of supercooled liquids above the glass 
transition, that the shear viscosity coefficient diverges at 
random close packing density, but with a larger exponent: 
rj ~ (n rcp — n) -/3 , (3 > 1. In our model of a crystallized 
shear flow we suggest a different approach. We will ac- 
commodate a recent finding of Luding et al. |l6j | that 
the shear viscosity coefficient diverges like (n* — n)^ 1 at 
a density 71* < n c : 



_ r l/2 



7? = 7?l 



(71* — n) d? 



(10) 



where 771 is a numerical factor of order unity, which is 
presently unknown. Divergence of the shear viscosity im- 
plies that the fluid is jammed on a macroscopic length 
scale. While a shear flow in a macroscopically jammed 
system is impossible, the system of hard spheres may still 
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have finite temperature, pressure, heat conduction, and 
collisional energy loss. 

There are arguments in the literature (see, e.g., San- 
tos et al. 01) that steady sheared states of granular 
fluids are intrinsically non-Newtonian (that is, require a 
description beyond Navier-Stokes order). Why are such 
effects not included in our description? The reason is 
that non-Newtonian effects arise when one treats the in- 
elasticity of particle collisions in a non-perturbative way 
(postulating that the Boltzmann equation, generalized to 
inelastic collisions, remains applicable for finite inelastic- 
ity). In this paper we work in the limit of nearly elastic 
particle collisions. In the first order in the inelasticity 
one does not need to take into account any inelasticity 
corrections to the transport coefficients (for example, to 
the shear viscosity). The only place where the inelastic- 
ity enters in this leading-order theory is the inelastic loss 
term T in the energy equation @. 



III. STEADY SHEAR FLOW CLOSE TO 
CRYSTALLIZATION 



addressed by Chou [Tg- He considered a model where 
inelastic spheres were driven by walls with attached half- 
spheres (bumpy walls.) Extending an earlier treatment 
by Richman [l!j, Chou calculated the heat flux and the 
slip velocity at the boundaries. He showed that the heat 
flux at the boundaries can be positive or negative, de- 
pending on whether the "slip work" is larger or smaller 
than the energy loss due to inelastic particle collisions 
with the walls. For some values of parameters the total 
heat flux to the boundaries vanishes • For simplicity, 
we will assume a vanishing heat flux and therefore pre- 
scribe dT/dy(y = 0) = dT/dy(y = H) = 0. The total 
number of particles is conserved, which yields a normal- 
ization condition for the density: n(y) dy — N, where 
N is the number of particles per unit area in the xz plane. 

Let us rescale the vertical coordinate y by the (a priori 
unknown) system height H, the horizontal velocity u by 
the upper plate velocity uq, the density n by the close- 
packing density n c , and the temperature T by the ratio 
Pq /n c of the constant applied pressure and close-packing 
density. Rewriting the equations in the scaled form, we 
obtain 



Now let us introduce the flow setting we will be dealing 
with in the rest of the paper. We consider a plane Cou- 
ette geometry. The model system is infinite in the hori- 
zontal (x) direction and driven by the upper wall y = H 
that moves horizontally with velocity uq . The lower wall 
y = is at rest (of course, in the absence of gravity the 
upper and lower walls are interchangeable). The height 
H of the layer in this setting is not fixed. Instead, the 
upper wall is maintained at a constant pressure Pq, like 
in a recent experiment by Gollub's group As already 
mentioned above, we assume a very dense flow, so that 
the volume fraction is close everywhere to that of close 
packing: (f> fcc - 4> < (f> fcc . 

As we consider a steady horizontal motion, the number 
density n, the granular temperature T and the horizon- 
tal velocity u depend only on the vertical coordinate y. 
Then it follows from the y-component of the momentum 
equation J3J) that the pressure is constant throughout the 
system: 



T 



p(y) = Po ■ 



(ii) 



Now we write down the n-component of the momentum 
equation J2J an d rewrite the energy balance equation J3J : 



d 
dy 



du 
dy 



= 0, 



d f dT 
dy \ d V 



du 



t?( — ) -T(n,T)=0, 



(12) 



where the constitutive relations are given by Eqs. 
(|10fl . Equations (JSJ-JT^J must be complemented by 
boundary conditions. We will assume rough walls and 
no-slip boundary conditions for the horizontal velocity: 
u{y — 0) = 0, u{y = H) = uq. The problem of eval- 
uation of the granular heat flux at the rough walls was 



Pi 



m MT 1 / 2 
K\ 1 — A — n 



d_ 

dy 



T l/2 



1-71 

du 



1 — A — n dy 

d_ f TV 2 dT 
dy \ 1 — n dy 

2 Ti 2i?T 3 / 2 

K\ 1 — 71 



1, 



= 0, 



(13) 



where M = n c u\jPQ is the square of the effective Mach 
number of the flow, and A = 1 — n* /n c is a small posi- 
tive numerical factor. For simplicity, the three presently 
unknown constants pi — 0(1), t] 1 /k 1 = 0(1) and 
Ti/ki = 0(1) are taken to be unity in the following. As 
the system height is unknown a priori, the scaled quan- 
tity R = (1 — r 2 ) H 2 /(2d 2 ) must be determined from 
the solution of the problem. From the first of Eqs. I|13|) 
we find a simple relation between the scaled density and 
temperature: T(y) = 1 — n[y). Then, introducing a con- 
venient auxiliary variable e(y) = [1 — rt(y)] 1 / 2 <C 1, we 
rewrite the remaining Eqs. (|13fl as 



(14) 




d e , , cA 

— + (c-R)e =0, 

dy e 



(15) 



where c is an unknown constant to be found from the so- 
lution of the problem. The boundary and normalization 
conditions are 

u(y = 0)=0, u(y = l) = l, 
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de , de , 

^^ = °) = ^^ = 1 ) = °' 



(16) 



where / = (1 — r 2 ) 1 / 2 iVeP/ 2. The steady flow equations 
O-tEJ include two (known) scaled parameters: M and 
/. There are five unknown parameters in the problem: 
c, R and three arbitrary constants which determine the 
solutions of the ordinary differential equations l|14|l and 
(|15|l for u(y) and e(y), respectively. Correspondingly, 
there are five conditions (|16|l to determine the unknown 
parameters. In the following, when solving the equations 
numerically, we put A = 0.05. 

Before we get into a detailed description of the steady 
flow solutions of Eqs. I)14 [l -I)16 [l . here is an overview of the 
results. Strikingly, there is an infinite number of steady 
flow solutions. They all can be parameterized by three 
parameters: the scaled numbers M and / and an integer 
number m = 1, 2, . . .. For a fixed m, there are three pos- 
sible types of solutions in different regions of the phase 
diagram (M, /), see Fig.|2] The simplest solution is the 
uniform, or linear shear flow that exists for any values of 
M and /. Here the velocity gradient, density and tem- 
perature are all constant. In region B (between the solid 
and dashed lines of Fig. EJ) , there is an additional one- 
phase solution: the one with nonlinear profiles of density, 
temperature and velocity. Finally, there is a multi-layer 
two-phase solution which exists in regions A and B, that 
is below the dashed line. The flow there is organized in 
several distinct layers. In the "solid-phase" layers the 
density is larger than the critical density n* at which the 
shear viscosity diverges. Therefore, these layers move as 
a whole with a velocity u = const, while the temperature 
and density profiles are non-trivial. In the "fluid phase" 
layers there is a mean flow with non-trivial profiles of 
the velocity, temperature and density. The density, tem- 
perature, heat flux and velocity are continuous at the 
interface between the layers. This means in particular, 
that there is no macroscopic flow in the bottom layer of 
the granulate, and the inelastic energy losses there are 
balanced by the conduction of heat from the (flowing) 
top layer. 

Now let us consider these solutions in more detail. We 
start with the simplest one: the uniform shear flow. Here 
e is independent of y, and Eq. ifPol) yields 



cA 

:-R' 

Substituting this value into Eq. Ijl4[) and integrating in 
y, we obtain a linear velocity profile 



u(y) = R , 



2A 



y- 



(17) 
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Using the normalization condition and the condition 
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FIG. 2: Hydro dynamic phase diagram of the steady flow so- 
lutions for m = 1 (the upper panel) and m = 2 (the lower 
panel) in terms of the scaled parameters / and M. In the 
regions A there are: a multi-layer flow and an (unstable) uni- 
form flow. In the regions B there are: a multi-layer flow, a 
non-uniform flow and an (unstable) uniform flow. In the re- 
gions C there is only a (stable) uniform flow. In the forbidden 
regions the granulate is not dense enough for our theory to be 
valid. Plausible bifurcation diagrams, corresponding to the 
dotted arrows / = 2, are shown, for the cases of m = 1 and 
m = 2 in Fig. 



u(y = 1) = 1, we obtain 
cA 



1 - 



f A P ■ 2XR2 

and c = R- 



c-R \AR M 
Solving these two algebraic equations, we calculate R 



(18) 



R 



2(1 -A) 2 ' 
M(l - A) 
P 



2M(1 - A)\ 1/2 



(19) 



and obtain the scaled velocity profile and the constant 
value of e: 



uo(y) 



v, 



A + 



(1 - A)fci 



1 + h + (1 + 2fci)V2 



1/2 



. (20) 
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where fci = M (1 — A) f~ 2 . Notice that e* always satis- 
fies the condition e > A 1 / 2 or, in the original variables, 
n < n*: the density of the uniform flow is smaller than 
the critical density at which the shear viscosity diverges. 
Note also that although the uniform shear flow solution 
[Eq. I|2()|l] formally exists everywhere, the assumption 
e<l demands fci <C 1. Therefore, the solution is valid 
for M <C / 2 / (l — A). This inequality breaks down in the 
forbidden regions of the phase diagram, see Fig. [5] 

The nonlinear solution, which exists in region B, can 
be found numerically. We used the following numerical 
procedure, realized in Matlab. Let us denote eo = e(y — 
0). For fixed R and /, we first solve Eq. (|T5)l by varying 
parameters eo and c and demanding de/dy(y = 1) = 
and the normalization condition. This procedure yields 
e(y). Then we solve Eq. I|14|) and find the velocity profile 
u(y) and the corresponding value of parameter M from 
the condition u(y = 1) = 1. 

An important insight into the nature of this solu- 
tion is provided by a mechanical analogy following from 
Eq. I|15|) . Let e be a "coordinate" , while the verti- 
cal coordinate y be "time". Then Eq. Q15[) describes 
a Newtonian particle oscillating in the potential well 
U(e) = (l/2)(c— R) e 2 — cAlne. The boundary conditions 
select a family of solutions. A typical solution includes 
an integer number to of halves of the oscillation period, 
and a completion of these oscillations takes the system 
a unit "time", mr/2 = 1 [sec Eq. <|21l) ]. where r is the 
oscillation period. Figure shows an example of particle 
in the potential well, while the resulting spatial profiles 
of e(y), u(y) and n(y) for a one-half of a full oscillation 
(to = 1), and for a full oscillation (to = 2), are shown in 
Fig.H 

At a fixed m, this solution exists in region B of the 
phase plane of parameters / and M. This region is lim- 
ited by two curves M = M(f) (see Fig. 0). The first 
curve is obtained as one approaches the bottom of the 
potential well. Here the oscillation amplitude (that is, 
the density contrast in the system) vanishes. Expand- 
ing the potential U(e) near the bottom of the potential 
well e = [cA/(c — -R)] 1 / 2 , we can calculate the oscillation 
period r: 



2 

TO 



R 



1/2 



(21) 



The bottom of the well corresponds to a uniform shear 
flow, but Eq. (|21|) must be obeyed arbitrarily close to the 
bottom of the well. Using Eq. (|2l")l and the second equa- 
tion of Eqs. HBJ), we obtain R = (tott/2) (Af/A) 1/2 and 
c = R + to 2 7t 2 /2. Substituting it into the first equation 
of Eqs. (|18|l . we obtain the lower boundary of region B 
in terms of / and M: 



f 



TO7T\ 



1/2 



1/4 



1 - A 



(AM) 



1/2 



(22) 



For m = 1 and to = 2 these curves are shown by the solid 
lines in Fig. [2 They determine the boundary between the 
regions A and B for each to. 
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FIG. 3: The mechanical analogy for Eq. flU} . Shown is a 
trajectory (the dashed line) of a classical particle in the po- 
tential well U (e) (the solid line, see text) which corresponds to 
regions B in Fig. El S ee Figure^Jfor the resulting spatial pro- 
files, which correspond to a one-half of the oscillation (m — 1) 
and a full oscillation (m = 2, here the potential changes as 
the parameter M is different.) The oscillating solutions exist 
only inside the potential well, and only for e(y — 0) > A 1 / 2 , 
when the density is everywhere below the critical density at 
which the shear viscosity diverges. The parameters for this 
figure are M = 0.5158 and / = 2. 




FIG. 4: The profiles of e(y) = [1 - n{y)] 1/2 (the upper row), 
the scaled velocity u(y) (the middle row), and the scaled den- 
sity n(y) (the bottom row), which correspond to region B 
in Figure |21 The profiles in the left column (respectively, in 
the right column) correspond to m — 1, a one-half of a full 
oscillation (respectively, to m = 2, a full oscillation) in the 
potential well, see Fig. The parameters are / = 2 and 
M = 0.5158 (the left column), and / = 2 and M = 0.1050 
(the right column). 



The second limiting curve (plotted by the dashed lines 
in Fig. |3 for to = 1 and to = 2) is obtained numerically 
by putting eo = A 1 / 2 (see Fig. |3J|- This curve determines 
the boundary between the regions B and C. 

Now let us consider the most interesting multi-layer 
two-phase family of solutions, parameterized by the same 
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number to > 1 that now can be associated with the num- 
ber of "solid" layers. Due to continuity of the velocity 
field, there is no macroscopic flow in the bottom layer 
of the granulate for any to, while the inelastic energy 
losses there are balanced by the heat conduction from 
the fluid-like top layer. A typical flow here consists of m 
zero-shear (solid-like) layers, each of which moving as a 
whole, separated by fluid-like regions. As we are deal- 
ing with a hard-sphere model, the granulate is fluidized 
within the "solid" layers: the granular temperature is 
non-zero, and there is energy flow through the bound- 
aries of the layers. Let us first consider the to — 1 so- 
lution. Here the only solid-like layer is at the bottom, 
and it is at rest. Putting c = in Eqs. (fl5|l and (fT4"jl we 
obtain e(y) = eo cosh(i? 1 / 2 ?/). The height of the bottom 
layer is determined from the condition e(h) — A 1 / 2 . At 
the interface between the two layers we demand continu- 
ity of the density (and, therefore, of the temperature), 
of the heat flux, and of the velocity. Similarly to the 
procedure performed in region B, we solve the problem 
numerically by shooting in two parameters eo and c (for 
the top layer) for fixed values of R and /, and then calcu- 
late the profiles and the respective value of the parameter 
M. Typical profiles e(y), u(y) and n(y) for m = 1 are 
shown in Fig. the left column. A multi-layer solution 
to > 1 is obtained in a similar way, by demanding that 
the velocity in any solid layer is constant, while the den- 
sity, temperature, heat flux and velocity are continuous 
at all interfaces between the layers. Typical profiles of 
e(y), u(y) and n(y) for to = 2 are shown in Fig. 03 right 
panel. 




FIG. 5: The profiles of e(y) (the upper row), the velocity 
u(y) (the middle row) and the density n(y) (the bottom row) 
for the two-phase solution, for m = 1 (the left column) and 
m = 2 (the right column). The solutions for the fluid (solid) 
layers are shown by the solid (dotted) lines. The parameters 
are / = 2 and M = 0.0492 (the left column), and / = 2 and 
M = 0.10148 (the right column.) 

At this point let us return, for a moment, to the case 
where the viscosity singularity is assumed to be at the 
same density n c as the rest of transport coefficient. In 



this case A = 0, and the only possible solution is a uni- 
form shear flow. Indeed, as e must be positive, the only 
acceptable solution of Eq. ((Tol with A = 0, which obeys 
the no-flux boundary conditions, is e = const and c = R. 
Then the normalization condition [the last of Eqs. I|16|l ] 
yields e 2 = 1 — f/R 1 ^ 2 . One can obtain this solution 
directly by putting A = in Eq. I|2(J|) . Therefore, the 
assumption of a nonzero A [the specific form of viscos- 
ity divergence, Eq. (|10|l ]. is crucial for the existence of 
non-trivial solutions. 

One can see that, at fixed to, two different kinds of 
steady flow solutions exist in region A, for the same val- 
ues of parameters M and /. Furthermore, three different 
kinds of steady flow solutions exist in region B, again 
for the same values of parameters M and / (see Fig. |2Jl . 
What is the selection rule for these solutions? We give a 
partial answer to this question in the next section by per- 
forming a linear stability analysis of the uniform dense 
shear flow. Then we suggest plausible bifurcation dia- 
grams of the system for different to. 



IV. LINEAR STABILITY AND BIFURCATIONS 

Our linear stability analysis of the uniform flow deals 
with the full set of (time-dependent) hydrodynamic equa- 
tions 0-© and constitutive relations (|5jl- (fTH . The 
details of linear stability analysis are shown in Ap- 
pendix. Adding a small perturbations to the uni- 
form shear flow and linearizing the equations, we fi- 
nally arrive at a quadratic characteristic equation for the 
growth/damping rate T as a function of parameters and 
the wave number k [Eq. (|A.4|) . see Appendix.] The two 
roots of this equation are real and correspond to two 
different collective modes of the system. One of them 
always decays. The other one is a purely growing mode 
for sufficiently small k, that is for long-wavelength per- 
turbations. At small k we obtain: 



r = 



2R 2 Xk 2 



M Le* (M + /ff/2) 



>0, 



(23) 



where L = (2M) 1/2 H/d. For sufficiently large k, T 
is negative: the heat conduction suppresses the insta- 
bility. The critical wave number for the instability is 
£;* = 2i?A 1 / 2 /M 1 / 2 . The dependence of the scaled growth 
rate T on the scaled wave number k is shown in Fig. [5] 
Using Eq. 119fl . we obtain 



f 2 A 1 / 2 



2(l-A) 2 M!/2 
M(l - A) 



f 2 



2M(1 - A) 
T 2 



1/2 



(24) 



Notice that the time scale separation, employed here for 
the reduction of the order of the dispersion equation (see 
Appendix), breaks down when e* approaches A 1 / 2 . As 
one can see from Eq. (|20|l . this happens, for a fixed 
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FIG. 6: The scaled growth rate F, calculated from Eq. 1A.4L 
versus the scaled wave number k. The long-wavelength modes 
with k < fc* (the value of k* is indicated by the asterisk) are 
unstable. The parameters / = 2, M = 0.2, A = 0.05, and 
L — 100 correspond to region A in Fig.|2K and to region C in 
Fig. UK 



/, when M becomes sufficiently small. Importantly, the 
final results {HJ| and {2D remain valid in the general 
case, as we obtained from the full, unreduced fourth- 
order dispersion equation for F(fc). 

The wave number k of the perturbation is quantized 
by the boundary conditions. Indeed, the velocity per- 
turbations must vanish at the upper and lower plates: 
m(y = 0, t) = m(y = = v 1 (y = 0,t) = v x {y = 
l,t) = 0. These conditions yield a discrete spectrum of 
wave numbers: k = nm, where to — 1,2,3,.... There- 
fore, the instability threshold is fc» = tt. This equation 
determines a curve f(M) on the (/, M) plane. One can 
see that this curve coincides with the curve f(M) deter- 
mined by Eq. (|22[1 for m = 1: the borderline between 
regions A and B in Fig. [21 Similarly, the threshold for 
to = 2 perturbation fc„ = 27r coincides with the curve 
f(M) determined by Eq. {22) for m = 2, and so on. This 
is not entirely surprising: as the instability is aperiodic, 
its threshold is provided by the marginal stability con- 
dition. It is convenient to characterize each steady flow 
solution by the maximum density contrast it predicts. In 
terms of the auxiliary variable e, we can define 5 as the 
maximum change in e throughout the system [for m = 1 
this gives S = e(y = 1) — e(y = 0).] We can calculate 
5 as a function of parameters / and M for each steady 
flow solution. Figure {7\ shows all the resulting branches 
of the solution (for m = 1 and m = 2) for a fixed /, while 
M serves as a control parameter. Furthermore, there is 
an infinite number of bifurcation diagrams for multi-layer 
solutions (for all integers m > 1), which can be computed 
in a similar way. 

In order to determine the bifurcation diagram of the 
system at a fixed m, one needs to perform a linear sta- 
bility analysis of each branch of the solution. Unfortu- 
nately, such an analysis is quite cumbersome for the so- 
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FIG. 7: A conjectured bifurcation diagrams of the steady 
flow for f — 2 and A = 0.05. Shown is the density contrast 
S versus M for different steady flow solutions for m = 1 (the 
upper panel) and m = 2 (the lower panel). The stable and 
unstable branches are indicated by the solid and dashed lines, 
respectively. We proved that the uniform flow (8 = 0) is 
unstable at M < M*. The dash-dotted line in the upper 
panel is the asymptote of the branch bifurcating from the 
point M = M»: S = 0.9856 . . . (M — M») 1/2 . 



lutions with 5 =/= 0. It is natural to assume that bifurca- 
tion, which occurs in the vicinity of the point M = M*, 
at which the uniform shear flow (the one with 5 = 0) 
changes from unstable [at M < M*(/)] to stable [at 
M > M*(/)], is an inverse pitchfork bifurcation. Based 
on this assumption, we obtain a simple bifurcation dia- 
gram, presented in Fig. [7] The stable solutions are in- 
dicated by solid lines, the unstable solutions by dashed 
lines. Let us consider the case of to = 1 and follow the 
uniform shear flow solution (for which 8 — 0) as M in- 
creases, starting from a small value, at fixed /. We pass 
from region A to region B of the hydrodynamic phase 
diagram (see the dotted arrow in Fig. [2J via an inverse 
pitchfork bifurcation. Here the uniform shear flow solu- 
tion becomes stable, while the bifurcating "second solu- 
tion" (the one with nonlinear density and velocity pro- 
files) must be unstable. Exploiting the mechanical anal- 
ogy (see Fig. 01, we can find the asymptote of the unsta- 
ble branch in the close vicinity of the bifurcation point: 
6 = A(f)(M - M*) 1 / 2 . Now let us move along the un- 
stable branch 6 ^ and increase M. At some critical 
value of M , which depends on /, we reach the border be- 
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tween regions B and C. Here the unstable branch ends, 
as for larger M (region C) the only possible solution is 
the uniform shear flow. There is, however, another solu- 
tion which exists at smaller M. Figure[7|shows this stable 
branch which corresponds to the two-phase solution in re- 
gions A and B (see Fig.|2Jl. This simple scenario predicts 
bistability at sufficiently small values of the parameter 
M, and a hysteresis on the interval [M*(/), Mi(/)], see 
Fig. Note that the parameter L does not affect the bi- 
furcation diagram, it only sets the time scales of transient 
motions in the system. 



V. SUMMARY AND DISCUSSION 

We have considered here some aspects of shear-induced 
crystallization in a dense but rapid mono-disperse gran- 
ular shear flow. We focused on a steady crystallized flow 
under a constant pressure and zero gravity. Assuming 
very high densities, n c — n <C n c , we employed a version 
of the Navier-Stokes hydrodynamics for inelastic hard 
sphere fluid with ad hoc constitutive relations based on 
the free volume argument. In contrast to earlier works 
on rapid granular shear flow, we assumed that the shear 
viscosity coefficient rj diverges at a density smaller than 
the close packing density n c , while the rest of the consti- 
tutive relations diverge at n = n c . We have determined 
the phase diagram of the steady flow in terms of three pa- 
rameters: the effective Mach number, the scaled inelastic 
energy loss parameter, and an integer number m. In a 
steady flow the viscous heating of the granulate is bal- 
anced by energy dissipation through inelastic collisions. 
This balance is achieved, in different parts of the phase 
diagram, in different ways, producing either a uniform 
shear flow (with constant velocity gradient, density and 
temperature), or a flow with nonlinear velocity, density 
and temperature profiles. In some regions in the phase 
diagram two or three different steady flow solutions are 
possible for the same values of the parameters. We per- 
formed a linear stability analysis of the uniform flow, and 
suggested a plausible bifurcation diagram of the flow at 
a fixed m, which predicts bi-stability and hysteresis. We 
are unable as yet to find a selection principle that would 
prefer certain steady state solutions out of a multitude of 
solutions at different m. This non-trivial selection prob- 
lem should be addressed in the future work. 

One of the predictions of this work is the existence 
(and, we conjecture, stability) of two-phase solutions. 
The simplest solution of this type consists of a zero-shear 
(solid-like) layer at the bottom and a flowing top layer. 
Though there is no mean flow in the bottom layer, the 
particles there undergo "thermal" motion, and the gran- 
ular temperature and pressure are non-zero. As a re- 
sult, there is energy transfer through the bottom layer. 
There also exist two-phase multi-layer solutions, where 
solid-like layers, each of which moving as a whole, are 
separated by fluid-like regions with nonlinear velocity, 
density and temperature profiles. The existence of these 



solutions is a direct consequence of our assumption that 
the coefficient of shear viscosity r\ in Eq. i|IU|) diverges at 
a density which is smaller than the close packing density 
n c . 

A comparison of our results with those of Alam and 
coworkers is in order now. Alam et al. |20| investigated 
a plane Couette shear flow of inelastic hard sphere fluid, 
in a system with a fixed height, in a wide interval of den- 
sities: from the dilute limit to the random close pack- 
ing density. They employed constitutive relations, all 
of which (including the viscosity) diverge at the random 
close packing density n r . They found instability of the 
uniform shear flow when the inelasticity of the particle 
collisions becomes large enough or, alternatively, when 
the (fixed) system height exceeds a critical value for a 
fixed inelasticity. The uniform shear flow instability con- 
sidered in our work is quite different: it requires viscosity 
divergence at a smaller density than the rest of constitu- 
tive relations. Indeed, Eqs. I|22[) and (|24|) show that this 
instability disappears when A goes to zero. Although the 
parameter A was identically zero in Ref. |20|. instability 
was observed. Where does the difference come from? To 
remind the reader, the constitutive relations that we used 
assume a very dense system. Respectively, our results are 
valid only at leading order in the parameter (n c — n)/n c . 
We checked that, if one takes into account only the lead- 
ing order terms in (n r — n)/n r in the equations of Alam 
et al. 20], the instability disappears. 

Its worth mentioning that the steady flow equations 
(|I4H - I|I6[) would not change if the crystal close packing 
density n c were replaced by the random close packing 
density n r . Such a formulation would follow from the 
assumption that the shear viscosity diverges at a density 
smaller than n r . It would predict a variety of steady state 
solutions on the metastable branch, in the vicinity of the 
random close packing. We hope that the basic assump- 
tions of our model (including the specific form of viscosity 
divergence) and its non-trivial predictions will be tested 
in MD simulations and in experiment on dense but rapid 
granular flow. Our work focused on a dense but rapid 
granular flow, assuming that the granulate is fully flu- 
idized. In experiment it should be easier to achieve this 
regime when the shear is very large, so that the effective 
Mach number is of order unity (see Fig. [7|) Most ex- 
periments with dense granular flows are performed with 
slow flows, where the effective Mach number is small. 
For example, in the experiment of Gollub's group 5] the 
parameter M was about 5 x 10~ 5 . In this regime the 
particles far from the moving boundary are in persistent 
contact with each other [5j, inter-particle friction is im- 
portant (see also [2l]]), and the model of inelastic hard 
spheres (and the Navier-Stokes hydrodynamics) is inap- 
plicable. 

In general, the model of inelastic hard sphere fluid is 
considered as a good approximation for dilute and mod- 
erately dense flows. Its validity range for dense flows is 
not well known ||. Indeed, the particle collision rate 
increases with the density, so the assumption of instan- 
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taneous collisions, intrinsic in the model of hard spheres, 
may become restrictive at high densities. Being aware 
of this limitation, we still believe that the model of in- 
elastic hard spheres can capture some of the physics of 
dense flow. Like in many other problems, it is useful to 
push this model to an extreme and analyze its predic- 
tions (some of which being quite unexpected as we have 
shown) in detail. We notice in passing that the two-phase 
flow predicted in this work resembles experimentally ob- 
served shear bands: localized regions of ordered granular 
flow, coexisting with essentially immobile solid-like re- 
gions 22]. 

The future experimental and theoretical work should 
elucidate the exact conditions under which shear-induced 
crystallization develops in granular flows. Though crys- 
tallization under shearing is well documented in MD sim- 
ulations 0,0, a recent experiment [23| showed that, un- 
der certain conditions, shearing can lead to disorder. 

Finally, we did not attempt to describe in this work the 
shear-induced crystallization process. Such a description 
is beyond the reach of theory as of present. A promising 
approach to this problem should deal, in addition to the 
hydrodynamic fields, with an order parameter field and 
its dynamics. For slow granular flows such a description 
is now emerging [24|. 
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APPENDIX: DETAILS OF LINEAR STABILITY 
ANALYSIS 

Let u and v be the velocity components in the x and 
y directions, respectively. We assume that the small per- 
turbations to the uniform steady shear flow do not break 
the symmetry of the flow in the x and z directions: 

n(y,t) = n + ni(y,t) , 
T(y,t)=T + T 1 (y,t), 
u(y,t) = u Q (y) + u 1 (y,t), 
v(y,t)=vi(y,t), (A.l) 

where index denotes the uniform steady shear flow, see 
Eqs. (|20(l . while Tq — 1 — no = e* 2 . Linearizing the hydro- 
dynamic equations with respect to small perturbations, 
we arrive at 



Ln 



n 



dvi 
dy 




vi 



dy 



Ln 



1/2 r a 
1-X-uq [dy 
dvi L dpi 
~dt = ~M ~dy 



2T 



"l 



1 - A 



n 



du Q 
dy 




(A.2) 



where L = (2 M) 1 ' 2 H/d, the time is measured in the 
units of H/uq, and the value of H is determined by the 
unperturbed flow, that is by the uniform steady flow so- 
lution. One can check that the new parameter L (which 
is absent in the steady state problem, but enters the lin- 
ear stability analysis), is fully determined by the scaled 
parameters / and M and by the restitution coefficient r. 

A further simplification employs time scale separation. 
Let us compare the characteristic time scales of the prob- 
lem. The acoustic time scale is t\ = H e* /(Po/^c) 1 ^ 2 , 
the heat conduction time scale is T2 = T\ (H/d), the en- 
ergy loss time scale is T3 = t\ [H/(dR)], and the vis- 
cous time is T4 = n (H/d) (e* 2 — A)/e* 2 . If the density 
is not too close to the density at which the shear vis- 
cosity diverges, that is if (H/d) (e* 2 — A)/e* 2 ^> 1, one 
can separate the different time scales and eliminate the 
(acoustic-like) fast modes. This is equivalent to the as- 
sumption that the perturbations evolve in pressure equi- 
librium with the surroundings |25j . Using the condition 
dpi/dt = instead of the full momentum equation [the 
third equation in Eqs. i|A.2|) ]. we obtain T\ — —n\. Then, 
differentiating the second equation of Eqs. I|A.2[) with re- 
spect to y and substituting dvi/dy from the first equation 
of Eqs. (|A.2|I . we arrive at 



dy 2 
dni 

M e* 
e 2- A 



dui 
A [ dy 
3(1 -e 2 ) 



dui 



ni 
"2e* 2 
1 



dy 



1 — (* 



2e« 



ni \ 


2 _ 


a; 


d 




dy 




3. 


+ 



1 dni 
* 2 dy 
Rni 



0, 



0. 



where we substituted duo/dy = 1 and To = 1 — no = e* 2 . 
Consider a single Fourier mode of the perturbation: 

ni (y, t) = h exp(rt + iky) , 
u i(y, t) = u exp(rt + iky) , 
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where k is the wave number. Looking for nontrivial so- where the coefficients are 
lutions, we obtain 

det(A)=0, (A.3) 

where the elements of the matrix A are T , 2 ft , 3 f 2 

a 2 = L k 1 



TLkf 2Re*k 3 



A 22 = 4 R k e* . 
Equation (| A. 3|> yields a quadratic equation for T: 

a 2 T 2 + ai r + a = 0, (A.4) 



2R 

M 2 M e* ' ai =4RLke* ( 1 



^ T 4i? 2 e*fc 2 i?fc 2 , , 



M 



3/ f ^* V M M/ 2 
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